#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Mar  1 12:28:48 2023

Reads in current and plasma data from 225 dayside "quiet" magnetopause crossings and 23 dayside EDR events.

Outputs event locations and current histograms comparing quiet and EDR current distributions. 

Also outputs a table of mean and median values for each current.

Note: "MPC" stands for magnetopause crossing in the code below, while "EDR" stands for electron diffusion region

@author: Jason Beedle
"""

import numpy as np
import matplotlib.pyplot as plt

#magnetopause crossing data storage class:
class magnetopauseCrossingData():
    
    def __init__(self,num):
        
        self.initial_num = num

def hist_reader(mpc,file):

  #Note, all currents are recoreded in spherical GSE coordiantes - see Beedle et. al 2022 for a full explanation
  Jr_curl   = [] ; Jp_curl   = [] ; Jt_curl   = []  #Curlometer current
  Jr_plasma = [] ; Jp_plasma = [] ; Jt_plasma = []  #FPI current
  Jr_para   = [] ; Jp_para   = [] ; Jt_para   = []  #Parallel Curlometer current
  Jr_perp   = [] ; Jp_perp   = [] ; Jt_perp   = []  #Perp Curlometer current
  Jr_totali = [] ; Jp_totali = [] ; Jt_totali = []  #Total ion diamagnetic current
  Jr_totale = [] ; Jp_totale = [] ; Jt_totale = []  #Total elec diamagnetic current
  
  Jr_gradNi = [] ; Jp_gradNi = [] ; Jt_gradNi = [] #Ion density gradient diamagnetic current component
  Jr_divTi  = [] ; Jp_divTi  = [] ; Jt_divTi  = [] #Ion temperature divergence diamagnetic current component
  
  Jr_gradNe = [] ; Jp_gradNe = [] ; Jt_gradNe = [] #Elec density gradient diamagnetic current component
  Jr_divTe  = [] ; Jp_divTe  = [] ; Jt_divTe  = [] #Elec temperature divergence diamagnetic current component
  
  Vi_r = [] ; Vi_p = [] ; Vi_t = [] #Ion Velocity
  Ve_r = [] ; Ve_p = [] ; Ve_t = [] #Electron Velocity
  
  B_x = [] ; B_y = [] ; B_z = [] #Magnetic Field - GSE coordinates
  
  Ti_para = [] ; Ti_perp = [] ; Ni = [] #Ion parallel and perpendicular temperature, and ion density
  
  #Opens data file
  with open (file, "r") as file:
        data = file.readlines()
        
        for line in data: 
            
            sections = line.split()
            
            #Reads out header line:
            if sections[0] == 'Jr_curl':
                print("first line")
            
            #Reads in and stores data:
            else:
                #Currents stored in micro A/m^2
                Jr_curl.append(float(sections[0]))
                Jp_curl.append(float(sections[1]))
                Jt_curl.append(float(sections[2]))
                
                Jr_para.append(float(sections[3]))
                Jp_para.append(float(sections[4]))
                Jt_para.append(float(sections[5]))
                
                Jr_perp.append(float(sections[6]))
                Jp_perp.append(float(sections[7]))
                Jt_perp.append(float(sections[8]))
                
                Jr_plasma.append(float(sections[9]))
                Jp_plasma.append(float(sections[10]))
                Jt_plasma.append(float(sections[11]))
            
                Jr_totali.append(float(sections[12]))
                Jp_totali.append(float(sections[13]))
                Jt_totali.append(float(sections[14]))
                
                Jr_gradNi.append(float(sections[15]))
                Jp_gradNi.append(float(sections[16]))
                Jt_gradNi.append(float(sections[17]))
                
                Jr_divTi.append(float(sections[18]))
                Jp_divTi.append(float(sections[19]))
                Jt_divTi.append(float(sections[20]))
                
                Jr_totale.append(float(sections[21]))
                Jp_totale.append(float(sections[22]))
                Jt_totale.append(float(sections[23]))
                
                Jr_gradNe.append(float(sections[24]))
                Jp_gradNe.append(float(sections[25]))
                Jt_gradNe.append(float(sections[26]))
                
                Jr_divTe.append(float(sections[27]))
                Jp_divTe.append(float(sections[28]))
                Jt_divTe.append(float(sections[29]))
                
                #km/s
                Vi_r.append(float(sections[30])) 
                Vi_p.append(float(sections[31]))
                Vi_t.append(float(sections[32]))
                
                Ve_r.append(float(sections[33]))
                Ve_p.append(float(sections[34]))
                Ve_t.append(float(sections[35]))
                
                #nT
                B_x.append(float(sections[36]))
                B_y.append(float(sections[37]))
                B_z.append(float(sections[38]))
                
                #eV
                Ti_para.append(float(sections[39]))
                Ti_perp.append(float(sections[40]))
                
                #cm^3
                Ni.append(float(sections[41]))
                
  #stores data in mpc class as arrays:
  mpc.Jr_curls = np.array(Jr_curl)
  mpc.Jp_curls = np.array(Jp_curl)
  mpc.Jt_curls = np.array(Jt_curl) 
  
  mpc.Jr_paras = np.array(Jr_para)
  mpc.Jp_paras = np.array(Jp_para)
  mpc.Jt_paras = np.array(Jt_para) 
  
  mpc.Jr_perps = np.array(Jr_perp)
  mpc.Jp_perps = np.array(Jp_perp)
  mpc.Jt_perps = np.array(Jt_perp) 
  
  mpc.Jr_plasmas = np.array(Jr_plasma)
  mpc.Jp_plasmas = np.array(Jp_plasma)
  mpc.Jt_plasmas = np.array(Jt_plasma) 
  
  mpc.Jr_totalis = np.array(Jr_totali)
  mpc.Jp_totalis = np.array(Jp_totali)
  mpc.Jt_totalis = np.array(Jt_totali)
  
  mpc.Jr_totales = np.array(Jr_totale)
  mpc.Jp_totales = np.array(Jp_totale)
  mpc.Jt_totales = np.array(Jt_totale)
  
  mpc.Vi_rs = np.array(Vi_r)
  mpc.Vi_ps = np.array(Vi_p)
  mpc.Vi_ts = np.array(Vi_t)
  
  mpc.Ve_rs = np.array(Ve_r)
  mpc.Ve_ps = np.array(Ve_p)
  mpc.Ve_ts = np.array(Ve_t)
  
  mpc.Jr_gradNis = np.array(Jr_gradNi)
  mpc.Jp_gradNis = np.array(Jp_gradNi)
  mpc.Jt_gradNis = np.array(Jt_gradNi)
  
  mpc.Jr_divTis = np.array(Jr_divTi)
  mpc.Jp_divTis = np.array(Jp_divTi)
  mpc.Jt_divTis = np.array(Jt_divTi)
  
  mpc.Jr_gradNes = np.array(Jr_gradNe)
  mpc.Jp_gradNes = np.array(Jp_gradNe)
  mpc.Jt_gradNes = np.array(Jt_gradNe)
  
  mpc.Jr_divTes = np.array(Jr_divTe)
  mpc.Jp_divTes = np.array(Jp_divTe)
  mpc.Jt_divTes = np.array(Jt_divTe)
  
  mpc.B_xs = np.array(B_x)
  mpc.B_ys = np.array(B_y)
  mpc.B_zs = np.array(B_z)
  
  mpc.Ti_paras = np.array(Ti_para)
  mpc.Ti_perps = np.array(Ti_perp)
  mpc.Nis = np.array(Ni)
  
  return  

def averaged_data_reader(mpc,file):

  #Initializes data arrays:
  MP_ID = [] ; MP_Start = [] ; MP_End = [] #MPC ID number and MPC start and end time
  
  MPC_date = [] #Year of magnetopause crossing
    
  CS_d = [] ; CS_Vn = []  #MPC duration (seconds) and normal velocity as determined by the deHoffman-Teller method
  
  Xmp = [] ; Ymp = [] ; Zmp = []  #Location of MPC based on Cartesian GSE coordinates (RE)
  thetaMP = [] ; phiMP = [] #Location of MPC based on our spherical GSE phi and theta angles (in degrees)
  
  CST_km = [] ; CST_di = [] ; CST_rgi = [] #MPC thickness in (km), scaled to di, and scaled to rgi 
  
  di = [] ; rgi = [] #Ion plasma scales

  #All current data is recorded in Spherical GSE coordiantes and consitutes averages of the current data points recorded over each event's MP crossing
  Jr_curl   = [] ; Jp_curl   = [] ; Jt_curl   = []  #Curlometer current
  Jr_para   = [] ; Jp_para   = [] ; Jt_para   = []  #Parallel Curlometer current
  Jr_plasma = [] ; Jp_plasma = [] ; Jt_plasma = []  #FPI current
  
  Jr_totali = [] ; Jp_totali = [] ; Jt_totali = []  #Total ion diamagnetic current
  Jr_divTi  = [] ; Jp_divTi  = [] ; Jt_divTi  = []  #Divergence T ion diamagnetic current component
  Jr_gradNi = [] ; Jp_gradNi = [] ; Jt_gradNi = []  #Gradient N ion diamagnetic current component
  
  Jr_totale = [] ; Jp_totale = [] ; Jt_totale = [] #Total elec diamagnetic current
  
  #Plasma data: averaged over all four MMS spacecraft during the magnetopause crossing
  Ni = [] ; Ti_para = [] ; Ti_perp = [] #Ion plasma parameters
  
  Vi_r = [] ; Vi_p = [] ; Vi_t = [] #Ion Velocity
  Ve_r = [] ; Ve_p = [] ; Ve_t = [] #Electron Velocity
  
  Bx = [] ; By = [] ; Bz = []  #Magnetic field in Cartesian GSE coordinates (nT)
  
  #Opens data file
  with open (file, "r") as file:
        data = file.readlines()
        
        for line in data: 
            
            sections = line.split()
            
            #Reads out header line:
            if sections[0] == 'MP_ID':
                print("first line")
            
            #Only reads in dayside values:
            elif float(sections[9]) > 50. or float(sections[9]) < -50.:
                print("flanks")
            
            #Reads in data to arrays
            else:
                MP_ID.append(sections[0])
                MP_Start.append(sections[1])
                MP_End.append(sections[2])
                
                MPC_date.append(sections[1][0:10])
                
                CS_d.append(float(sections[3]))
                CS_Vn.append(float(sections[4]))
            
                Xmp.append(float(sections[5]))
                Ymp.append(float(sections[6]))
                Zmp.append(float(sections[7]))
            
                thetaMP.append(float(sections[8]))
                phiMP.append(float(sections[9]))
            
                CST_km.append(float(sections[10]))
                CST_di.append(float(sections[11]))
                CST_rgi.append(float(sections[12]))
                
                di.append(float(sections[13]))
                rgi.append(float(sections[14]))
                
                Jr_curl.append(float(sections[15]))
                Jp_curl.append(float(sections[16]))
                Jt_curl.append(float(sections[17]))
            
                Jr_para.append(float(sections[18]))
                Jp_para.append(float(sections[19]))
                Jt_para.append(float(sections[20]))
                
                Jr_plasma.append(float(sections[21]))
                Jp_plasma.append(float(sections[22]))
                Jt_plasma.append(float(sections[23]))
            
                Jr_totali.append(float(sections[24]))
                Jp_totali.append(float(sections[25]))
                Jt_totali.append(float(sections[26]))
            
                Jr_gradNi.append(float(sections[27]))
                Jp_gradNi.append(float(sections[28]))
                Jt_gradNi.append(float(sections[29]))
            
                Jr_divTi.append(float(sections[30]))
                Jp_divTi.append(float(sections[31]))
                Jt_divTi.append(float(sections[32]))
            
                Jr_totale.append(float(sections[33]))
                Jp_totale.append(float(sections[34]))
                Jt_totale.append(float(sections[35]))
            
                Ti_para.append(float(sections[36]))
                Ti_perp.append(float(sections[37]))
                Ni.append(float(sections[38]))
            
                Vi_r.append(float(sections[39]))
                Vi_p.append(float(sections[40]))
                Vi_t.append(float(sections[41]))
                
                Ve_r.append(float(sections[42]))
                Ve_p.append(float(sections[43]))
                Ve_t.append(float(sections[44]))
                
                Bx.append(float(sections[45]))
                By.append(float(sections[46]))
                Bz.append(float(sections[47]))
                
  #stores data in mpc class:
  mpc.MP_ID = np.array(MP_ID) 
  
  mpc.num_events = len(MP_ID)
  
  mpc.MPC_date = np.array(MPC_date)

  mpc.MP_Start = MP_Start
  mpc.MP_End   = MP_End           
                
  mpc.CS_d  = CS_d
  mpc.CS_Vn = CS_Vn
                
  mpc.Xmp = Xmp
  mpc.Ymp = Ymp
  mpc.Zmp = Zmp
  
  mpc.thetaMP = thetaMP
  mpc.phiMP   = phiMP
  
  mpc.CST_km  = np.array(CST_km)
  mpc.CST_di  = CST_di
  mpc.CST_rgi = CST_rgi
  
  mpc.di = di
  mpc.rgi = rgi
   
  mpc.Jr_curl = np.array(Jr_curl)
  mpc.Jp_curl = np.array(Jp_curl)
  mpc.Jt_curl = np.array(Jt_curl) 
  
  mpc.Jr_para = np.array(Jr_para)
  mpc.Jp_para = np.array(Jp_para)
  mpc.Jt_para = np.array(Jt_para) 
  
  mpc.Jr_plasma = np.array(Jr_plasma)
  mpc.Jp_plasma = np.array(Jp_plasma)
  mpc.Jt_plasma = np.array(Jt_plasma) 
  
  mpc.Jr_totali = np.array(Jr_totali)
  mpc.Jp_totali = np.array(Jp_totali)
  mpc.Jt_totali = np.array(Jt_totali)
  
  mpc.Jr_gradNi = np.array(Jr_gradNi)
  mpc.Jp_gradNi = np.array(Jp_gradNi)
  mpc.Jt_gradNi = np.array(Jt_gradNi)
  
  mpc.Jr_divTi = np.array(Jr_divTi)
  mpc.Jp_divTi = np.array(Jp_divTi)
  mpc.Jt_divTi = np.array(Jt_divTi)
  
  mpc.Jr_totale = np.array(Jr_totale)
  mpc.Jp_totale = np.array(Jp_totale)
  mpc.Jt_totale = np.array(Jt_totale)
  
  mpc.Ti_para = np.array(Ti_para)
  mpc.Ti_perp = np.array(Ti_perp)
  
  mpc.Ni = np.array(Ni) 
  
  mpc.Vi_r = np.array(Vi_r)
  mpc.Vi_p = np.array(Vi_p)
  mpc.Vi_t = np.array(Vi_t)
  
  mpc.Ve_r = np.array(Ve_r)
  mpc.Ve_p = np.array(Ve_p)
  mpc.Ve_t = np.array(Ve_t)
  
  mpc.Bx = Bx 
  mpc.By = By
  mpc.Bz = Bz
  
  return  

#Plots 6 panel normalized histograms - used for ion and electron velocity data
def histPlotter6Panel(one_edr,two_edr,three_edr,four_edr,five_edr,six_edr,one_mpc,two_mpc,three_mpc,four_mpc,five_mpc,six_mpc,titles,fileTitle,xtitle,one_b,two_b,three_b,bin_size):
    
    figHist,axHist = plt.subplots(2,3,figsize=(30,20))
    
    fontsizes = 35
    
    axHist[0,0].set_title(titles[0],fontsize=fontsizes)
    axHist[0,1].set_title(titles[1],fontsize=fontsizes)
    axHist[0,2].set_title(titles[2],fontsize=fontsizes)
    
    axHist[1,0].set_title(titles[3],fontsize=fontsizes)
    axHist[1,1].set_title(titles[4],fontsize=fontsizes)
    axHist[1,2].set_title(titles[5],fontsize=fontsizes)
    
    #Creates uniform bins to use for all six histograms
    bins = np.histogram(np.hstack((one_edr[0],two_edr[0],three_edr[0],four_edr[0],five_edr[0],six_edr[0],one_mpc[0],two_mpc[0],three_mpc[0],four_mpc[0],five_mpc[0],six_mpc[0])), bins=bin_size)[1]
    
    #Creates distributions for EDR data - in blue
    axHist[0,0].hist(one_edr[0],bins,density=True,color='blue',label='N: %i'%(len(one_edr[0])) +'\n'+'Mean: %0.2f'%(one_edr[1]) +'\n'+'Median: %0.1f \n'%(one_edr[2]))
    axHist[0,1].hist(two_edr[0],bins,density=True,color='blue',label='N: %i'%(len(two_edr[0])) +'\n'+'Mean: %0.2f'%(two_edr[1]) +'\n'+'Median: %0.1f \n'%(two_edr[2]))
    axHist[0,2].hist(three_edr[0],bins,density=True,color='blue',label='N: %i'%(len(three_edr[0])) +'\n'+'Mean: %0.2f'%(three_edr[1]) +'\n'+'Median: %0.1f \n'%(three_edr[2]))
    
    #Creates distributions for quiet data - in red
    axHist[0,0].hist(one_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(one_mpc[0])) +'\n'+'Mean: %0.2f'%(one_mpc[1]) +'\n'+'Median: %0.1f '%(one_mpc[2]))
    axHist[0,1].hist(two_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(two_mpc[0])) +'\n'+'Mean: %0.2f'%(two_mpc[1]) +'\n'+'Median: %0.1f '%(two_mpc[2]))
    axHist[0,2].hist(three_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(three_mpc[0])) +'\n'+'Mean: %0.2f'%(three_mpc[1]) +'\n'+'Median: %0.1f '%(three_mpc[2]))
    
    #Creates distributions for EDR data - in blue
    axHist[1,0].hist(four_edr[0],bins,density=True,color='blue',label='N: %i'%(len(four_edr[0])) +'\n'+'Mean: %0.2f'%(four_edr[1]) +'\n'+'Median: %0.1f \n'%(four_edr[2]))
    axHist[1,1].hist(five_edr[0],bins,density=True,color='blue',label='N: %i'%(len(five_edr[0])) +'\n'+'Mean: %0.2f'%(five_edr[1]) +'\n'+'Median: %0.1f \n'%(five_edr[2]))
    axHist[1,2].hist(six_edr[0],bins,density=True,color='blue',label='N: %i'%(len(six_edr[0])) +'\n'+'Mean: %0.2f'%(six_edr[1]) +'\n'+'Median: %0.1f \n'%(six_edr[2]))
   
    #Creates distributions for quiet data - in red
    axHist[1,0].hist(four_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(four_mpc[0])) +'\n'+'Mean: %0.2f'%(four_mpc[1]) +'\n'+'Median: %0.1f '%(four_mpc[2]))
    axHist[1,1].hist(five_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(five_mpc[0])) +'\n'+'Mean: %0.2f'%(five_mpc[1]) +'\n'+'Median: %0.1f '%(five_mpc[2]))
    axHist[1,2].hist(six_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(six_mpc[0])) +'\n'+'Mean: %0.2f'%(six_mpc[1]) +'\n'+'Median: %0.1f '%(six_mpc[2]))
    
    axHist[0,0].set_xlim(one_b)
    axHist[0,1].set_xlim(two_b)
    axHist[0,2].set_xlim(three_b)
    
    axHist[1,0].set_xlim(one_b)
    axHist[1,1].set_xlim(two_b)
    axHist[1,2].set_xlim(three_b)
    
    axHist[0,0].legend(fontsize=(fontsizes/2),loc='upper right')
    axHist[0,1].legend(fontsize=(fontsizes/2),loc='upper left')
    axHist[0,2].legend(fontsize=(fontsizes/2),loc='upper right')
    
    axHist[1,0].legend(fontsize=(fontsizes/2),loc='upper right')
    axHist[1,1].legend(fontsize=(fontsizes/2),loc='upper left')
    axHist[1,2].legend(fontsize=(fontsizes/2),loc='upper right')
    
    axHist[1,1].set_xlabel(xtitle,fontsize=fontsizes)
    
    axHist[0,0].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[0,1].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[0,2].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    
    axHist[1,0].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[1,1].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[1,2].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    
    axHist[0,0].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[0,1].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[0,2].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    
    axHist[1,0].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[1,1].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[1,2].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    
    axHist[0,0].get_shared_y_axes().join(axHist[0,0],axHist[0,1],axHist[0,2],axHist[1,0],axHist[1,1],axHist[1,2])
    
    figHist.savefig(fileTitle+'_Hist'+'.png')
    
    return

#Plots 9 panel normalized histograms - used for curlometer, ion, and electron current data
def histPlotter9Panel(one_edr,two_edr,three_edr,four_edr,five_edr,six_edr,seven_edr,eight_edr,nine_edr,one_mpc,two_mpc,three_mpc,four_mpc,five_mpc,six_mpc,seven_mpc,eight_mpc,nine_mpc,titles,fileTitle,xtitle,one_b,two_b,three_b,bin_size):
    
    figHist,axHist = plt.subplots(3,3,figsize=(30,30))
    
    fontsizes = 35
    
    axHist[0,0].set_title(titles[0],fontsize=fontsizes)
    axHist[0,1].set_title(titles[1],fontsize=fontsizes)
    axHist[0,2].set_title(titles[2],fontsize=fontsizes)
    
    axHist[1,0].set_title(titles[3],fontsize=fontsizes)
    axHist[1,1].set_title(titles[4],fontsize=fontsizes)
    axHist[1,2].set_title(titles[5],fontsize=fontsizes)
    
    axHist[2,0].set_title(titles[6],fontsize=fontsizes)
    axHist[2,1].set_title(titles[7],fontsize=fontsizes)
    axHist[2,2].set_title(titles[8],fontsize=fontsizes)
    
    #Creates uniform bins to use for all nine histograms
    bins = np.histogram(np.hstack((one_edr[0],two_edr[0],three_edr[0],four_edr[0],five_edr[0],six_edr[0],seven_edr[0],eight_edr[0],nine_edr[0],one_mpc[0],two_mpc[0],three_mpc[0],four_mpc[0],five_mpc[0],six_mpc[0],seven_mpc[0],eight_mpc[0],nine_mpc[0])), bins=bin_size)[1]
    
    #Creates distributions for EDR data - in blue
    axHist[0,0].hist(one_edr[0],bins,density=True,color='blue',label='N: %i'%(len(one_edr[0])) +'\n'+'Mean: %0.1f'%(one_edr[1]) +'\n'+'Median: %0.1f \n'%(one_edr[2]))
    axHist[0,1].hist(two_edr[0],bins,density=True,color='blue',label='N: %i'%(len(two_edr[0])) +'\n'+'Mean: %0.1f'%(two_edr[1]) +'\n'+'Median: %0.1f \n'%(two_edr[2]))
    axHist[0,2].hist(three_edr[0],bins,density=True,color='blue',label='N: %i'%(len(three_edr[0])) +'\n'+'Mean: %0.1f'%(three_edr[1]) +'\n'+'Median: %0.1f \n'%(three_edr[2]))
    
    #Creates distributions for quiet data - in red
    axHist[0,0].hist(one_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(one_mpc[0])) +'\n'+'Mean: %0.1f'%(one_mpc[1]) +'\n'+'Median: %0.1f '%(one_mpc[2]))
    axHist[0,1].hist(two_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(two_mpc[0])) +'\n'+'Mean: %0.1f'%(two_mpc[1]) +'\n'+'Median: %0.1f '%(two_mpc[2]))
    axHist[0,2].hist(three_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(three_mpc[0])) +'\n'+'Mean: %0.1f'%(three_mpc[1]) +'\n'+'Median: %0.1f'%(three_mpc[2]))
    
    #Creates distributions for EDR data - in blue
    axHist[1,0].hist(four_edr[0],bins,density=True,color='blue',label='N: %i'%(len(four_edr[0])) +'\n'+'Mean: %0.1f'%(four_edr[1]) +'\n'+'Median: %0.1f \n'%(four_edr[2]))
    axHist[1,1].hist(five_edr[0],bins,density=True,color='blue',label='N: %i'%(len(five_edr[0])) +'\n'+'Mean: %0.1f'%(five_edr[1]) +'\n'+'Median: %0.1f \n'%(five_edr[2]))
    axHist[1,2].hist(six_edr[0],bins,density=True,color='blue',label='N: %i'%(len(six_edr[0])) +'\n'+'Mean: %0.1f'%(six_edr[1]) +'\n'+'Median: %0.1f \n'%(six_edr[2]))
   
    #Creates distributions for quiet data - in red
    axHist[1,0].hist(four_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(four_mpc[0])) +'\n'+'Mean: %0.1f'%(four_mpc[1]) +'\n'+'Median: %0.1f '%(four_mpc[2]))
    axHist[1,1].hist(five_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(five_mpc[0])) +'\n'+'Mean: %0.1f'%(five_mpc[1]) +'\n'+'Median: %0.1f '%(five_mpc[2]))
    axHist[1,2].hist(six_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(six_mpc[0])) +'\n'+'Mean: %0.1f'%(six_mpc[1]) +'\n'+'Median: %0.1f '%(six_mpc[2]))
    
    #Creates distributions for EDR data - in blue
    axHist[2,0].hist(seven_edr[0],bins,density=True,color='blue',label='N: %i'%(len(seven_edr[0])) +'\n'+'Mean: %0.1f'%(seven_edr[1]) +'\n'+'Median: %0.1f \n'%(seven_edr[2]))
    axHist[2,1].hist(eight_edr[0],bins,density=True,color='blue',label='N: %i'%(len(eight_edr[0])) +'\n'+'Mean: %0.1f'%(eight_edr[1]) +'\n'+'Median: %0.1f \n'%(eight_edr[2]))
    axHist[2,2].hist(nine_edr[0],bins,density=True,color='blue',label='N: %i'%(len(nine_edr[0])) +'\n'+'Mean: %0.1f'%(nine_edr[1]) +'\n'+'Median: %0.1f \n'%(nine_edr[2]))
   
    #Creates distributions for quiet data - in red
    axHist[2,0].hist(seven_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(seven_mpc[0])) +'\n'+'Mean: %0.1f'%(seven_mpc[1]) +'\n'+'Median: %0.1f '%(seven_mpc[2]))
    axHist[2,1].hist(eight_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(eight_mpc[0])) +'\n'+'Mean: %0.1f'%(eight_mpc[1]) +'\n'+'Median: %0.1f '%(eight_mpc[2]))
    axHist[2,2].hist(nine_mpc[0],bins=bins,density=True,color='red',alpha = 0.4,label='N: %i'%(len(nine_mpc[0])) +'\n'+'Mean: %0.1f'%(nine_mpc[1]) +'\n'+'Median: %0.1f '%(nine_mpc[2]))
    
    
    axHist[0,0].set_xlim(one_b)
    axHist[0,1].set_xlim(two_b)
    axHist[0,2].set_xlim(three_b)
    
    axHist[1,0].set_xlim(one_b)
    axHist[1,1].set_xlim(two_b)
    axHist[1,2].set_xlim(three_b)
    
    axHist[2,0].set_xlim(one_b)
    axHist[2,1].set_xlim(two_b)
    axHist[2,2].set_xlim(three_b)
    
    axHist[0,0].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[0,1].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[0,2].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    
    axHist[1,0].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[1,1].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[1,2].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    
    axHist[2,0].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[2,1].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    axHist[2,2].legend(fontsize=(fontsizes/2)+5,loc='upper right')
    
    axHist[0,0].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[0,1].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[0,2].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    
    axHist[1,0].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[1,1].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[1,2].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    
    axHist[2,0].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[2,1].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    axHist[2,2].tick_params(axis='x', labelsize=(fontsizes/2)+5)
    
    axHist[0,0].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[0,1].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[0,2].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    
    axHist[1,0].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[1,1].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[1,2].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    
    axHist[2,0].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[2,1].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    axHist[2,2].tick_params(axis='y', labelsize=(fontsizes/2)+5)
    
    axHist[0,0].get_shared_y_axes().join(axHist[0,0],axHist[0,1],axHist[0,2],axHist[1,0],axHist[1,1],axHist[1,2],axHist[2,0],axHist[2,1],axHist[2,2])
    
    axHist[2,1].set_xlabel(xtitle,fontsize=fontsizes)
    
    figHist.savefig(fileTitle+'_Hist'+'.png')
    
    return

#Takes in EDR and Quiet Hist and gets their each distributions mean and median values, whose data is then passed to the plotting functions
def createHistograms(edr_hist,quiet_hist):
    
    #EDR data:
        
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_curls_edr = edr_hist.Jr_curls*1000.
    Jp_curls_edr = edr_hist.Jp_curls*1000.
    Jt_curls_edr = edr_hist.Jt_curls*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_curls_edr_mean   = np.nanmean(Jr_curls_edr)
    Jr_curls_edr_median = np.nanmedian(Jr_curls_edr)
    Jr_curls_edr_std    = np.nanstd(Jr_curls_edr)
    
    Jp_curls_edr_mean   = np.nanmean(Jp_curls_edr)
    Jp_curls_edr_median = np.nanmedian(Jp_curls_edr)
    Jp_curls_edr_std    = np.nanstd(Jp_curls_edr)
    
    Jt_curls_edr_mean   = np.nanmean(Jt_curls_edr)
    Jt_curls_edr_median = np.nanmedian(Jt_curls_edr)
    Jt_curls_edr_std    = np.nanstd(Jt_curls_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_plasmas_edr = edr_hist.Jr_plasmas[~np.isnan(edr_hist.Jr_plasmas)]*1000.
    Jp_plasmas_edr = edr_hist.Jp_plasmas[~np.isnan(edr_hist.Jp_plasmas)]*1000.
    Jt_plasmas_edr = edr_hist.Jt_plasmas[~np.isnan(edr_hist.Jt_plasmas)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_plasmas_edr_mean   = np.nanmean(Jr_plasmas_edr)
    Jr_plasmas_edr_median = np.nanmedian(Jr_plasmas_edr)
    Jr_plasmas_edr_std    = np.nanstd(Jr_plasmas_edr)
    
    Jp_plasmas_edr_mean   = np.nanmean(Jp_plasmas_edr)
    Jp_plasmas_edr_median = np.nanmedian(Jp_plasmas_edr)
    Jp_plasmas_edr_std    = np.nanstd(Jp_plasmas_edr)
    
    Jt_plasmas_edr_mean   = np.nanmean(Jt_plasmas_edr)
    Jt_plasmas_edr_median = np.nanmedian(Jt_plasmas_edr)
    Jt_plasmas_edr_std    = np.nanstd(Jt_plasmas_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_paras_edr = edr_hist.Jr_paras*1000.
    Jp_paras_edr = edr_hist.Jp_paras*1000.
    Jt_paras_edr = edr_hist.Jt_paras*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_paras_edr_mean   = np.nanmean(Jr_paras_edr)
    Jr_paras_edr_median = np.nanmedian(Jr_paras_edr)
    Jr_paras_edr_std    = np.nanstd(Jr_paras_edr)
    
    Jp_paras_edr_mean   = np.nanmean(Jp_paras_edr)
    Jp_paras_edr_median = np.nanmedian(Jp_paras_edr)
    Jp_paras_edr_std    = np.nanstd(Jp_paras_edr)
    
    Jt_paras_edr_mean   = np.nanmean(Jt_paras_edr)
    Jt_paras_edr_median = np.nanmedian(Jt_paras_edr)
    Jt_paras_edr_std    = np.nanstd(Jt_paras_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_perps_edr = edr_hist.Jr_perps*1000.
    Jp_perps_edr = edr_hist.Jp_perps*1000.
    Jt_perps_edr = edr_hist.Jt_perps*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_perps_edr_mean   = np.nanmean(Jr_perps_edr)
    Jr_perps_edr_median = np.nanmedian(Jr_perps_edr)
    Jr_perps_edr_std    = np.nanstd(Jr_perps_edr)
    
    Jp_perps_edr_mean   = np.nanmean(Jp_perps_edr)
    Jp_perps_edr_median = np.nanmedian(Jp_perps_edr)
    Jp_perps_edr_std    = np.nanstd(Jp_perps_edr)
    
    Jt_perps_edr_mean   = np.nanmean(Jt_perps_edr)
    Jt_perps_edr_median = np.nanmedian(Jt_perps_edr)
    Jt_perps_edr_std    = np.nanstd(Jt_perps_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totalis_edr = edr_hist.Jr_totalis[~np.isnan(edr_hist.Jr_totalis)]*1000.
    Jp_totalis_edr = edr_hist.Jp_totalis[~np.isnan(edr_hist.Jp_totalis)]*1000.
    Jt_totalis_edr = edr_hist.Jt_totalis[~np.isnan(edr_hist.Jt_totalis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_totalis_edr_mean   = np.nanmean(Jr_totalis_edr)
    Jr_totalis_edr_median = np.nanmedian(Jr_totalis_edr)
    Jr_totalis_edr_std    = np.nanstd(Jr_totalis_edr)
    
    Jp_totalis_edr_mean   = np.nanmean(Jp_totalis_edr)
    Jp_totalis_edr_median = np.nanmedian(Jp_totalis_edr)
    Jp_totalis_edr_std    = np.nanstd(Jp_totalis_edr)
    
    Jt_totalis_edr_mean   = np.nanmean(Jt_totalis_edr)
    Jt_totalis_edr_median = np.nanmedian(Jt_totalis_edr)
    Jt_totalis_edr_std    = np.nanstd(Jt_totalis_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totales_edr = edr_hist.Jr_totales*1000.
    Jp_totales_edr = edr_hist.Jp_totales*1000.
    Jt_totales_edr = edr_hist.Jt_totales*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_totales_edr_mean   = np.nanmean(Jr_totales_edr)
    Jr_totales_edr_median = np.nanmedian(Jr_totales_edr)
    Jr_totales_edr_std    = np.nanstd(Jr_totales_edr)
    
    Jp_totales_edr_mean   = np.nanmean(Jp_totales_edr)
    Jp_totales_edr_median = np.nanmedian(Jp_totales_edr)
    Jp_totales_edr_std    = np.nanstd(Jp_totales_edr)
    
    Jt_totales_edr_mean   = np.nanmean(Jt_totales_edr)
    Jt_totales_edr_median = np.nanmedian(Jt_totales_edr)
    Jt_totales_edr_std    = np.nanstd(Jt_totales_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNis_edr = edr_hist.Jr_gradNis[~np.isnan(edr_hist.Jr_gradNis)]*1000.
    Jp_gradNis_edr = edr_hist.Jp_gradNis[~np.isnan(edr_hist.Jp_gradNis)]*1000.
    Jt_gradNis_edr = edr_hist.Jt_gradNis[~np.isnan(edr_hist.Jt_gradNis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_gradNis_edr_mean   = np.nanmean(Jr_gradNis_edr)
    Jr_gradNis_edr_median = np.nanmedian(Jr_gradNis_edr)
    Jr_gradNis_edr_std    = np.nanstd(Jr_gradNis_edr)
    
    Jp_gradNis_edr_mean   = np.nanmean(Jp_gradNis_edr)
    Jp_gradNis_edr_median = np.nanmedian(Jp_gradNis_edr)
    Jp_gradNis_edr_std    = np.nanstd(Jp_gradNis_edr)
    
    Jt_gradNis_edr_mean   = np.nanmean(Jt_gradNis_edr)
    Jt_gradNis_edr_median = np.nanmedian(Jt_gradNis_edr)
    Jt_gradNis_edr_std    = np.nanstd(Jt_gradNis_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTis_edr = edr_hist.Jr_divTis[~np.isnan(edr_hist.Jr_divTis)]*1000.
    Jp_divTis_edr = edr_hist.Jp_divTis[~np.isnan(edr_hist.Jp_divTis)]*1000.
    Jt_divTis_edr = edr_hist.Jt_divTis[~np.isnan(edr_hist.Jt_divTis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_divTis_edr_mean   = np.nanmean(Jr_divTis_edr)
    Jr_divTis_edr_median = np.nanmedian(Jr_divTis_edr)
    Jr_divTis_edr_std    = np.nanstd(Jr_divTis_edr)
    
    Jp_divTis_edr_mean   = np.nanmean(Jp_divTis_edr)
    Jp_divTis_edr_median = np.nanmedian(Jp_divTis_edr)
    Jp_divTis_edr_std    = np.nanstd(Jp_divTis_edr)
    
    Jt_divTis_edr_mean   = np.nanmean(Jt_divTis_edr)
    Jt_divTis_edr_median = np.nanmedian(Jt_divTis_edr)
    Jt_divTis_edr_std    = np.nanstd(Jt_divTis_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNes_edr = edr_hist.Jr_gradNes[~np.isnan(edr_hist.Jr_gradNes)]*1000.
    Jp_gradNes_edr = edr_hist.Jp_gradNes[~np.isnan(edr_hist.Jp_gradNes)]*1000.
    Jt_gradNes_edr = edr_hist.Jt_gradNes[~np.isnan(edr_hist.Jt_gradNes)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_gradNes_edr_mean   = np.nanmean(Jr_gradNes_edr)
    Jr_gradNes_edr_median = np.nanmedian(Jr_gradNes_edr)
    Jr_gradNes_edr_std    = np.nanstd(Jr_gradNes_edr)
    
    Jp_gradNes_edr_mean   = np.nanmean(Jp_gradNes_edr)
    Jp_gradNes_edr_median = np.nanmedian(Jp_gradNes_edr)
    Jp_gradNes_edr_std    = np.nanstd(Jp_gradNes_edr)
    
    Jt_gradNes_edr_mean   = np.nanmean(Jt_gradNes_edr)
    Jt_gradNes_edr_median = np.nanmedian(Jt_gradNes_edr)
    Jt_gradNes_edr_std    = np.nanstd(Jt_gradNes_edr)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTes_edr = edr_hist.Jr_divTes[~np.isnan(edr_hist.Jr_divTes)]*1000.
    Jp_divTes_edr = edr_hist.Jp_divTes[~np.isnan(edr_hist.Jp_divTes)]*1000.
    Jt_divTes_edr = edr_hist.Jt_divTes[~np.isnan(edr_hist.Jt_divTes)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_divTes_edr_mean   = np.nanmean(Jr_divTes_edr)
    Jr_divTes_edr_median = np.nanmedian(Jr_divTes_edr)
    Jr_divTes_edr_std    = np.nanstd(Jr_divTes_edr)
    
    Jp_divTes_edr_mean   = np.nanmean(Jp_divTes_edr)
    Jp_divTes_edr_median = np.nanmedian(Jp_divTes_edr)
    Jp_divTes_edr_std    = np.nanstd(Jp_divTes_edr)
    
    Jt_divTes_edr_mean   = np.nanmean(Jt_divTes_edr)
    Jt_divTes_edr_median = np.nanmedian(Jt_divTes_edr)
    Jt_divTes_edr_std    = np.nanstd(Jt_divTes_edr)
    
    #Initializes ion velocities - reads out nan values
    Vi_r_edr = edr_hist.Vi_rs[~np.isnan(edr_hist.Vi_rs)]
    Vi_p_edr = edr_hist.Vi_ps[~np.isnan(edr_hist.Vi_ps)]
    Vi_t_edr = edr_hist.Vi_ts[~np.isnan(edr_hist.Vi_ts)]
    
    #Calculates the mean, median, and standard deviation of each velocity component
    Vi_r_edr_mean   = np.nanmean(Vi_r_edr)
    Vi_r_edr_median = np.nanmedian(Vi_r_edr)
    Vi_r_edr_std    = np.nanstd(Vi_r_edr)
    
    Vi_p_edr_mean   = np.nanmean(Vi_p_edr)
    Vi_p_edr_median = np.nanmedian(Vi_p_edr)
    Vi_p_edr_std    = np.nanstd(Vi_p_edr)
    
    Vi_t_edr_mean   = np.nanmean(Vi_t_edr)
    Vi_t_edr_median = np.nanmedian(Vi_t_edr)
    Vi_t_edr_std    = np.nanstd(Vi_t_edr)
    
    #Initializes elec velocities - reads out nan values
    Ve_r_edr = edr_hist.Ve_rs[~np.isnan(edr_hist.Ve_rs)]
    Ve_p_edr = edr_hist.Ve_ps[~np.isnan(edr_hist.Ve_ps)]
    Ve_t_edr = edr_hist.Ve_ts[~np.isnan(edr_hist.Ve_ts)]
    
    #Calculates the mean, median, and standard deviation of each velocity component
    Ve_r_edr_mean   = np.nanmean(Ve_r_edr)
    Ve_r_edr_median = np.nanmedian(Ve_r_edr)
    Ve_r_edr_std    = np.nanstd(Ve_r_edr)
    
    Ve_p_edr_mean   = np.nanmean(Ve_p_edr)
    Ve_p_edr_median = np.nanmedian(Ve_p_edr)
    Ve_p_edr_std    = np.nanstd(Ve_p_edr)
    
    Ve_t_edr_mean   = np.nanmean(Ve_t_edr)
    Ve_t_edr_median = np.nanmedian(Ve_t_edr)
    Ve_t_edr_std    = np.nanstd(Ve_t_edr)    
    
    #Quiet data:
        
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_curls_quiet = quiet_hist.Jr_curls*1000.
    Jp_curls_quiet = quiet_hist.Jp_curls*1000.
    Jt_curls_quiet = quiet_hist.Jt_curls*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_curls_quiet_mean   = np.nanmean(Jr_curls_quiet)
    Jr_curls_quiet_median = np.nanmedian(Jr_curls_quiet)
    Jr_curls_quiet_std    = np.nanstd(Jr_curls_quiet)
    
    Jp_curls_quiet_mean   = np.nanmean(Jp_curls_quiet)
    Jp_curls_quiet_median = np.nanmedian(Jp_curls_quiet)
    Jp_curls_quiet_std    = np.nanstd(Jp_curls_quiet)
    
    Jt_curls_quiet_mean   = np.nanmean(Jt_curls_quiet)
    Jt_curls_quiet_median = np.nanmedian(Jt_curls_quiet)
    Jt_curls_quiet_std    = np.nanstd(Jt_curls_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_plasmas_quiet = quiet_hist.Jr_plasmas[~np.isnan(quiet_hist.Jr_plasmas)]*1000.
    Jp_plasmas_quiet = quiet_hist.Jp_plasmas[~np.isnan(quiet_hist.Jp_plasmas)]*1000.
    Jt_plasmas_quiet = quiet_hist.Jt_plasmas[~np.isnan(quiet_hist.Jt_plasmas)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_plasmas_quiet_mean   = np.nanmean(Jr_plasmas_quiet)
    Jr_plasmas_quiet_median = np.nanmedian(Jr_plasmas_quiet)
    Jr_plasmas_quiet_std    = np.nanstd(Jr_plasmas_quiet)
    
    Jp_plasmas_quiet_mean   = np.nanmean(Jp_plasmas_quiet)
    Jp_plasmas_quiet_median = np.nanmedian(Jp_plasmas_quiet)
    Jp_plasmas_quiet_std    = np.nanstd(Jp_plasmas_quiet)
    
    Jt_plasmas_quiet_mean   = np.nanmean(Jt_plasmas_quiet)
    Jt_plasmas_quiet_median = np.nanmedian(Jt_plasmas_quiet)
    Jt_plasmas_quiet_std    = np.nanstd(Jt_plasmas_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_paras_quiet = quiet_hist.Jr_paras*1000.
    Jp_paras_quiet = quiet_hist.Jp_paras*1000.
    Jt_paras_quiet = quiet_hist.Jt_paras*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_paras_quiet_mean   = np.nanmean(Jr_paras_quiet)
    Jr_paras_quiet_median = np.nanmedian(Jr_paras_quiet)
    Jr_paras_quiet_std    = np.nanstd(Jr_paras_quiet)
    
    Jp_paras_quiet_mean   = np.nanmean(Jp_paras_quiet)
    Jp_paras_quiet_median = np.nanmedian(Jp_paras_quiet)
    Jp_paras_quiet_std    = np.nanstd(Jp_paras_quiet)
    
    Jt_paras_quiet_mean   = np.nanmean(Jt_paras_quiet)
    Jt_paras_quiet_median = np.nanmedian(Jt_paras_quiet)
    Jt_paras_quiet_std    = np.nanstd(Jt_paras_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_perps_quiet = quiet_hist.Jr_perps*1000.
    Jp_perps_quiet = quiet_hist.Jp_perps*1000.
    Jt_perps_quiet = quiet_hist.Jt_perps*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_perps_quiet_mean   = np.nanmean(Jr_perps_quiet)
    Jr_perps_quiet_median = np.nanmedian(Jr_perps_quiet)
    Jr_perps_quiet_std    = np.nanstd(Jr_perps_quiet)
    
    Jp_perps_quiet_mean   = np.nanmean(Jp_perps_quiet)
    Jp_perps_quiet_median = np.nanmedian(Jp_perps_quiet)
    Jp_perps_quiet_std    = np.nanstd(Jp_perps_quiet)
    
    Jt_perps_quiet_mean   = np.nanmean(Jt_perps_quiet)
    Jt_perps_quiet_median = np.nanmedian(Jt_perps_quiet)
    Jt_perps_quiet_std    = np.nanstd(Jt_perps_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totalis_quiet = quiet_hist.Jr_totalis[~np.isnan(quiet_hist.Jr_totalis)]*1000.
    Jp_totalis_quiet = quiet_hist.Jp_totalis[~np.isnan(quiet_hist.Jp_totalis)]*1000.
    Jt_totalis_quiet = quiet_hist.Jt_totalis[~np.isnan(quiet_hist.Jt_totalis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_totalis_quiet_mean   = np.nanmean(Jr_totalis_quiet)
    Jr_totalis_quiet_median = np.nanmedian(Jr_totalis_quiet)
    Jr_totalis_quiet_std    = np.nanstd(Jr_totalis_quiet)
    
    Jp_totalis_quiet_mean   = np.nanmean(Jp_totalis_quiet)
    Jp_totalis_quiet_median = np.nanmedian(Jp_totalis_quiet)
    Jp_totalis_quiet_std    = np.nanstd(Jp_totalis_quiet)
    
    Jt_totalis_quiet_mean   = np.nanmean(Jt_totalis_quiet)
    Jt_totalis_quiet_median = np.nanmedian(Jt_totalis_quiet)
    Jt_totalis_quiet_std    = np.nanstd(Jt_totalis_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totales_quiet = quiet_hist.Jr_totales*1000.
    Jp_totales_quiet = quiet_hist.Jp_totales*1000.
    Jt_totales_quiet = quiet_hist.Jt_totales*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_totales_quiet_mean   = np.nanmean(Jr_totales_quiet)
    Jr_totales_quiet_median = np.nanmedian(Jr_totales_quiet)
    Jr_totales_quiet_std    = np.nanstd(Jr_totales_quiet)
    
    Jp_totales_quiet_mean   = np.nanmean(Jp_totales_quiet)
    Jp_totales_quiet_median = np.nanmedian(Jp_totales_quiet)
    Jp_totales_quiet_std    = np.nanstd(Jp_totales_quiet)
    
    Jt_totales_quiet_mean   = np.nanmean(Jt_totales_quiet)
    Jt_totales_quiet_median = np.nanmedian(Jt_totales_quiet)
    Jt_totales_quiet_std    = np.nanstd(Jt_totales_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNis_quiet = quiet_hist.Jr_gradNis[~np.isnan(quiet_hist.Jr_gradNis)]*1000.
    Jp_gradNis_quiet = quiet_hist.Jp_gradNis[~np.isnan(quiet_hist.Jp_gradNis)]*1000.
    Jt_gradNis_quiet = quiet_hist.Jt_gradNis[~np.isnan(quiet_hist.Jt_gradNis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_gradNis_quiet_mean   = np.nanmean(Jr_gradNis_quiet)
    Jr_gradNis_quiet_median = np.nanmedian(Jr_gradNis_quiet)
    Jr_gradNis_quiet_std    = np.nanstd(Jr_gradNis_quiet)
    
    Jp_gradNis_quiet_mean   = np.nanmean(Jp_gradNis_quiet)
    Jp_gradNis_quiet_median = np.nanmedian(Jp_gradNis_quiet)
    Jp_gradNis_quiet_std    = np.nanstd(Jp_gradNis_quiet)
    
    Jt_gradNis_quiet_mean   = np.nanmean(Jt_gradNis_quiet)
    Jt_gradNis_quiet_median = np.nanmedian(Jt_gradNis_quiet)
    Jt_gradNis_quiet_std    = np.nanstd(Jt_gradNis_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTis_quiet = quiet_hist.Jr_divTis[~np.isnan(quiet_hist.Jr_divTis)]*1000.
    Jp_divTis_quiet = quiet_hist.Jp_divTis[~np.isnan(quiet_hist.Jp_divTis)]*1000.
    Jt_divTis_quiet = quiet_hist.Jt_divTis[~np.isnan(quiet_hist.Jt_divTis)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_divTis_quiet_mean   = np.nanmean(Jr_divTis_quiet)
    Jr_divTis_quiet_median = np.nanmedian(Jr_divTis_quiet)
    Jr_divTis_quiet_std    = np.nanstd(Jr_divTis_quiet)
    
    Jp_divTis_quiet_mean   = np.nanmean(Jp_divTis_quiet)
    Jp_divTis_quiet_median = np.nanmedian(Jp_divTis_quiet)
    Jp_divTis_quiet_std    = np.nanstd(Jp_divTis_quiet)
    
    Jt_divTis_quiet_mean   = np.nanmean(Jt_divTis_quiet)
    Jt_divTis_quiet_median = np.nanmedian(Jt_divTis_quiet)
    Jt_divTis_quiet_std    = np.nanstd(Jt_divTis_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNes_quiet = quiet_hist.Jr_gradNes[~np.isnan(quiet_hist.Jr_gradNes)]*1000.
    Jp_gradNes_quiet = quiet_hist.Jp_gradNes[~np.isnan(quiet_hist.Jp_gradNes)]*1000.
    Jt_gradNes_quiet = quiet_hist.Jt_gradNes[~np.isnan(quiet_hist.Jt_gradNes)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_gradNes_quiet_mean   = np.nanmean(Jr_gradNes_quiet)
    Jr_gradNes_quiet_median = np.nanmedian(Jr_gradNes_quiet)
    Jr_gradNes_quiet_std    = np.nanstd(Jr_gradNes_quiet)
    
    Jp_gradNes_quiet_mean   = np.nanmean(Jp_gradNes_quiet)
    Jp_gradNes_quiet_median = np.nanmedian(Jp_gradNes_quiet)
    Jp_gradNes_quiet_std    = np.nanstd(Jp_gradNes_quiet)
    
    Jt_gradNes_quiet_mean   = np.nanmean(Jt_gradNes_quiet)
    Jt_gradNes_quiet_median = np.nanmedian(Jt_gradNes_quiet)
    Jt_gradNes_quiet_std    = np.nanstd(Jt_gradNes_quiet)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTes_quiet = quiet_hist.Jr_divTes[~np.isnan(quiet_hist.Jr_divTes)]*1000.
    Jp_divTes_quiet = quiet_hist.Jp_divTes[~np.isnan(quiet_hist.Jp_divTes)]*1000.
    Jt_divTes_quiet = quiet_hist.Jt_divTes[~np.isnan(quiet_hist.Jt_divTes)]*1000.
    
    #Calculates the mean, median, and standard deviation of each current component
    Jr_divTes_quiet_mean   = np.nanmean(Jr_divTes_quiet)
    Jr_divTes_quiet_median = np.nanmedian(Jr_divTes_quiet)
    Jr_divTes_quiet_std    = np.nanstd(Jr_divTes_quiet)
    
    Jp_divTes_quiet_mean   = np.nanmean(Jp_divTes_quiet)
    Jp_divTes_quiet_median = np.nanmedian(Jp_divTes_quiet)
    Jp_divTes_quiet_std    = np.nanstd(Jp_divTes_quiet)
    
    Jt_divTes_quiet_mean   = np.nanmean(Jt_divTes_quiet)
    Jt_divTes_quiet_median = np.nanmedian(Jt_divTes_quiet)
    Jt_divTes_quiet_std    = np.nanstd(Jt_divTes_quiet)
    
    #Initializes ion velocities - reads out nan values
    Vi_r_quiet = quiet_hist.Vi_rs[~np.isnan(quiet_hist.Vi_rs)]
    Vi_p_quiet = quiet_hist.Vi_ps[~np.isnan(quiet_hist.Vi_ps)]
    Vi_t_quiet = quiet_hist.Vi_ts[~np.isnan(quiet_hist.Vi_ts)]
    
    #Calculates the mean, median, and standard deviation of each velocity component
    Vi_r_quiet_mean   = np.nanmean(Vi_r_quiet)
    Vi_r_quiet_median = np.nanmedian(Vi_r_quiet)
    Vi_r_quiet_std    = np.nanstd(Vi_r_quiet)
    
    Vi_p_quiet_mean   = np.nanmean(Vi_p_quiet)
    Vi_p_quiet_median = np.nanmedian(Vi_p_quiet)
    Vi_p_quiet_std    = np.nanstd(Vi_p_quiet)
    
    Vi_t_quiet_mean   = np.nanmean(Vi_t_quiet)
    Vi_t_quiet_median = np.nanmedian(Vi_t_quiet)
    Vi_t_quiet_std    = np.nanstd(Vi_t_quiet)
    
    #Initializes elec velocities - reads out nan values
    Ve_r_quiet = quiet_hist.Ve_rs[~np.isnan(quiet_hist.Ve_rs)]
    Ve_p_quiet = quiet_hist.Ve_ps[~np.isnan(quiet_hist.Ve_ps)]
    Ve_t_quiet = quiet_hist.Ve_ts[~np.isnan(quiet_hist.Ve_ts)]
    
    #Calculates the mean, median, and standard deviation of each velocity component
    Ve_r_quiet_mean   = np.nanmean(Ve_r_quiet)
    Ve_r_quiet_median = np.nanmedian(Ve_r_quiet)
    Ve_r_quiet_std    = np.nanstd(Ve_r_quiet)
    
    Ve_p_quiet_mean   = np.nanmean(Ve_p_quiet)
    Ve_p_quiet_median = np.nanmedian(Ve_p_quiet)
    Ve_p_quiet_std    = np.nanstd(Ve_p_quiet)
    
    Ve_t_quiet_mean   = np.nanmean(Ve_t_quiet)
    Ve_t_quiet_median = np.nanmedian(Ve_t_quiet)
    Ve_t_quiet_std    = np.nanstd(Ve_t_quiet)    
    
    #Plots Velocity Histograms
    histPlotter6Panel([Vi_r_edr,Vi_r_edr_mean,Vi_r_edr_median],[Vi_p_edr,Vi_p_edr_mean,Vi_p_edr_median],[Vi_t_edr,Vi_t_edr_mean,Vi_t_edr_median],[Ve_r_edr,Ve_r_edr_mean,Ve_r_edr_median],[Ve_p_edr,Ve_p_edr_mean,Ve_p_edr_median],[Ve_t_edr,Ve_t_edr_mean,Ve_t_edr_median],[Vi_r_quiet,Vi_r_quiet_mean,Vi_r_quiet_median],[Vi_p_quiet,Vi_p_quiet_mean,Vi_p_quiet_median],[Vi_t_quiet,Vi_t_quiet_mean,Vi_t_quiet_median],[Ve_r_quiet,Ve_r_quiet_mean,Ve_r_quiet_median],[Ve_p_quiet,Ve_p_quiet_mean,Ve_p_quiet_median],[Ve_t_quiet,Ve_t_quiet_mean,Ve_t_quiet_median],[r'$V_{iR}$',r'$V_{i \phi}$',r'$V_{i \theta}$',r'$V_{eR}$',r'$V_{e \phi}$',r'$V_{e \theta}$'],"Velocities","Velocity (km/s)",[-500,500],[-1100,600],[-750,1000],250)
    
    #Plots Ion Diamagnetic Current Histograms
    histPlotter9Panel([Jr_totalis_edr,Jr_totalis_edr_mean,Jr_totalis_edr_median],[Jp_totalis_edr,Jp_totalis_edr_mean,Jp_totalis_edr_median],[Jt_totalis_edr,Jt_totalis_edr_mean,Jt_totalis_edr_median],[Jr_gradNis_edr,Jr_gradNis_edr_mean,Jr_gradNis_edr_median],[Jp_gradNis_edr,Jp_gradNis_edr_mean,Jp_gradNis_edr_median],[Jt_gradNis_edr,Jt_gradNis_edr_mean,Jt_gradNis_edr_median],[Jr_divTis_edr,Jr_divTis_edr_mean,Jr_divTis_edr_median],[Jp_divTis_edr,Jp_divTis_edr_mean,Jp_divTis_edr_median],[Jt_divTis_edr,Jt_divTis_edr_mean,Jt_divTis_edr_median],[Jr_totalis_quiet,Jr_totalis_quiet_mean,Jr_totalis_quiet_median],[Jp_totalis_quiet,Jp_totalis_quiet_mean,Jp_totalis_quiet_median],[Jt_totalis_quiet,Jt_totalis_quiet_mean,Jt_totalis_quiet_median],[Jr_gradNis_quiet,Jr_gradNis_quiet_mean,Jr_gradNis_quiet_median],[Jp_gradNis_quiet,Jp_gradNis_quiet_mean,Jp_gradNis_quiet_median],[Jt_gradNis_quiet,Jt_gradNis_quiet_mean,Jt_gradNis_quiet_median],[Jr_divTis_quiet,Jr_divTis_quiet_mean,Jr_divTis_quiet_median],[Jp_divTis_quiet,Jp_divTis_quiet_mean,Jp_divTis_quiet_median],[Jt_divTis_quiet,Jt_divTis_quiet_mean,Jt_divTis_quiet_median],[r'$J_{R \ dia \ Total_i}$',r'$J_{\phi \ dia \ Total_i}$',r'$J_{\theta \ dia \ Total_i}$',r'$J_{R \ dia \ \nabla N_i}$',r'$J_{\phi \ dia \ \nabla N_i}$',r'$J_{\theta \ dia \ \nabla N_i}$',r'$J_{R \ dia \ \nabla \cdot T_i}$',r'$J_{\phi \ dia \ \nabla \cdot T_i}$',r'$J_{\theta \ dia \ \nabla \cdot T_i}$'],"Jion_dia_Currents","Current Density "+r'$ (nA/m^2)$',[-900,1500],[-900,1500],[-900,1500],250)
    
    #Plots Electron Diamagnetic Current Histograms
    histPlotter9Panel([Jr_totales_edr,Jr_totales_edr_mean,Jr_totales_edr_median],[Jp_totales_edr,Jp_totales_edr_mean,Jp_totales_edr_median],[Jt_totales_edr,Jt_totales_edr_mean,Jt_totales_edr_median],[Jr_gradNes_edr,Jr_gradNes_edr_mean,Jr_gradNes_edr_median],[Jp_gradNes_edr,Jp_gradNes_edr_mean,Jp_gradNes_edr_median],[Jt_gradNes_edr,Jt_gradNes_edr_mean,Jt_gradNes_edr_median],[Jr_divTes_edr,Jr_divTes_edr_mean,Jr_divTes_edr_median],[Jp_divTes_edr,Jp_divTes_edr_mean,Jp_divTes_edr_median],[Jt_divTes_edr,Jt_divTes_edr_mean,Jt_divTes_edr_median],[Jr_totales_quiet,Jr_totales_quiet_mean,Jr_totales_quiet_median],[Jp_totales_quiet,Jp_totales_quiet_mean,Jp_totales_quiet_median],[Jt_totales_quiet,Jt_totales_quiet_mean,Jt_totales_quiet_median],[Jr_gradNes_quiet,Jr_gradNes_quiet_mean,Jr_gradNes_quiet_median],[Jp_gradNes_quiet,Jp_gradNes_quiet_mean,Jp_gradNes_quiet_median],[Jt_gradNes_quiet,Jt_gradNes_quiet_mean,Jt_gradNes_quiet_median],[Jr_divTes_quiet,Jr_divTes_quiet_mean,Jr_divTes_quiet_median],[Jp_divTes_quiet,Jp_divTes_quiet_mean,Jp_divTes_quiet_median],[Jt_divTes_quiet,Jt_divTes_quiet_mean,Jt_divTes_quiet_median],[r'$J_{R \ dia \ Total_e}$',r'$J_{\phi \ dia \ Total_e}$',r'$J_{\theta \ dia \ Total_e}$',r'$J_{R \ dia \ \nabla N_e}$',r'$J_{\phi \ dia \ \nabla N_e}$',r'$J_{\theta \ dia \ \nabla N_e}$',r'$J_{R \ dia \ \nabla \cdot T_e}$',r'$J_{\phi \ dia \ \nabla \cdot T_e}$',r'$J_{\theta \ dia \ \nabla \cdot T_e}$'],"Jelec_dia_Currents","Current Density "+r'$ (nA/m^2)$',[-250,400],[-250,400],[-250,400],250)
    
    #Plots Curlometer Current Histograms
    histPlotter9Panel([Jr_curls_edr,Jr_curls_edr_mean,Jr_curls_edr_median],[Jp_curls_edr,Jp_curls_edr_mean,Jp_curls_edr_median],[Jt_curls_edr,Jt_curls_edr_mean,Jt_curls_edr_median],[Jr_perps_edr,Jr_perps_edr_mean,Jr_perps_edr_median],[Jp_perps_edr,Jp_perps_edr_mean,Jp_perps_edr_median],[Jt_perps_edr,Jt_perps_edr_mean,Jt_perps_edr_median],[Jr_paras_edr,Jr_paras_edr_mean,Jr_paras_edr_median],[Jp_paras_edr,Jp_paras_edr_mean,Jp_paras_edr_median],[Jt_paras_edr,Jt_paras_edr_mean,Jt_paras_edr_median],[Jr_curls_quiet,Jr_curls_quiet_mean,Jr_curls_quiet_median],[Jp_curls_quiet,Jp_curls_quiet_mean,Jp_curls_quiet_median],[Jt_curls_quiet,Jt_curls_quiet_mean,Jt_curls_quiet_median],[Jr_perps_quiet,Jr_perps_quiet_mean,Jr_perps_quiet_median],[Jp_perps_quiet,Jp_perps_quiet_mean,Jp_perps_quiet_median],[Jt_perps_quiet,Jt_perps_quiet_mean,Jt_perps_quiet_median],[Jr_paras_quiet,Jr_paras_quiet_mean,Jr_paras_quiet_median],[Jp_paras_quiet,Jp_paras_quiet_mean,Jp_paras_quiet_median],[Jt_paras_quiet,Jt_paras_quiet_mean,Jt_paras_quiet_median],[r'$J_{R \ curl}$',r'$J_{\phi \ curl}$',r'$J_{\theta \ curl}$',r'$J_{R \ curl_\perp}$',r'$J_{\phi \ curl_\perp}$',r'$J_{\theta \ curl_\perp}$',r'$J_{R \ curl_\parallel}$',r'$J_{\phi \ curl_\parallel}$',r'$J_{\theta \ curl_\parallel}$'],"Jcurl_Currents","Current Density "+r'$ (nA/m^2)$',[-900,1500],[-900,1500],[-900,1500],250)
    
    return

#Calcualtes the standard error of a given array using the standard deviation
def standErr(arr):
        
        num = len(arr)
        
        std = np.nanstd(arr)
        
        stanError = std/(np.sqrt(num))
            
        return stanError    

#calculates the mean, median, and standard error values of each current and outputs them to a table matching Table 1 in Beedle et al. 2023
def data_table(mpc,hist):
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_curl = hist.Jr_curls*1000.
    Jp_curl = hist.Jp_curls*1000.
    Jt_curl = hist.Jt_curls*1000.
    
    #Calculates the mean value
    Jr_curl_mean = np.nanmean(Jr_curl)
    Jp_curl_mean = np.nanmean(Jp_curl)
    Jt_curl_mean = np.nanmean(Jt_curl)
    
    #Calculates the median value
    Jr_curl_median = np.nanmedian(Jr_curl)
    Jp_curl_median = np.nanmedian(Jp_curl)
    Jt_curl_median = np.nanmedian(Jt_curl)
    
    #Calculates the standard error value
    Jr_curl_err = standErr(Jr_curl)
    Jp_curl_err = standErr(Jp_curl)
    Jt_curl_err = standErr(Jt_curl)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_perp = hist.Jr_perps*1000.
    Jp_perp = hist.Jp_perps*1000.
    Jt_perp = hist.Jt_perps*1000.
    
    Jr_perp_mean = np.nanmean(Jr_perp)
    Jp_perp_mean = np.nanmean(Jp_perp)
    Jt_perp_mean = np.nanmean(Jt_perp)
    
    Jr_perp_median = np.nanmedian(Jr_perp)
    Jp_perp_median = np.nanmedian(Jp_perp)
    Jt_perp_median = np.nanmedian(Jt_perp)
    
    Jr_perp_err = standErr(Jr_perp)
    Jp_perp_err = standErr(Jp_perp)
    Jt_perp_err = standErr(Jt_perp)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_para = hist.Jr_paras*1000.
    Jp_para = hist.Jp_paras*1000.
    Jt_para = hist.Jt_paras*1000.
    
    Jr_para_mean = np.nanmean(Jr_para)
    Jp_para_mean = np.nanmean(Jp_para)
    Jt_para_mean = np.nanmean(Jt_para)
    
    Jr_para_median = np.nanmedian(Jr_para)
    Jp_para_median = np.nanmedian(Jp_para)
    Jt_para_median = np.nanmedian(Jt_para)
    
    Jr_para_err = standErr(Jr_para)
    Jp_para_err = standErr(Jp_para)
    Jt_para_err = standErr(Jt_para)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totali = hist.Jr_totalis*1000.
    Jp_totali = hist.Jp_totalis*1000.
    Jt_totali = hist.Jt_totalis*1000.
    
    Jr_totali_mean = np.nanmean(Jr_totali)
    Jp_totali_mean = np.nanmean(Jp_totali)
    Jt_totali_mean = np.nanmean(Jt_totali)
    
    Jr_totali_median = np.nanmedian(Jr_totali)
    Jp_totali_median = np.nanmedian(Jp_totali)
    Jt_totali_median = np.nanmedian(Jt_totali)
    
    Jr_totali_err = standErr(Jr_totali)
    Jp_totali_err = standErr(Jp_totali)
    Jt_totali_err = standErr(Jt_totali)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNi = hist.Jr_gradNis*1000.
    Jp_gradNi = hist.Jp_gradNis*1000.
    Jt_gradNi = hist.Jt_gradNis*1000.
    
    Jr_gradNi_mean = np.nanmean(Jr_gradNi)
    Jp_gradNi_mean = np.nanmean(Jp_gradNi)
    Jt_gradNi_mean = np.nanmean(Jt_gradNi)
    
    Jr_gradNi_median = np.nanmedian(Jr_gradNi)
    Jp_gradNi_median = np.nanmedian(Jp_gradNi)
    Jt_gradNi_median = np.nanmedian(Jt_gradNi)
    
    Jr_gradNi_err = standErr(Jr_gradNi)
    Jp_gradNi_err = standErr(Jp_gradNi)
    Jt_gradNi_err = standErr(Jt_gradNi)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTi = hist.Jr_divTis*1000.
    Jp_divTi = hist.Jp_divTis*1000.
    Jt_divTi = hist.Jt_divTis*1000.
    
    Jr_divTi_mean = np.nanmean(Jr_divTi)
    Jp_divTi_mean = np.nanmean(Jp_divTi)
    Jt_divTi_mean = np.nanmean(Jt_divTi)
    
    Jr_divTi_median = np.nanmedian(Jr_divTi)
    Jp_divTi_median = np.nanmedian(Jp_divTi)
    Jt_divTi_median = np.nanmedian(Jt_divTi)
    
    Jr_divTi_err = standErr(Jr_divTi)
    Jp_divTi_err = standErr(Jp_divTi)
    Jt_divTi_err = standErr(Jt_divTi)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_totale = hist.Jr_totales*1000.
    Jp_totale = hist.Jp_totales*1000.
    Jt_totale = hist.Jt_totales*1000.
    
    Jr_totale_mean = np.nanmean(Jr_totale)
    Jp_totale_mean = np.nanmean(Jp_totale)
    Jt_totale_mean = np.nanmean(Jt_totale)
    
    Jr_totale_median = np.nanmedian(Jr_totale)
    Jp_totale_median = np.nanmedian(Jp_totale)
    Jt_totale_median = np.nanmedian(Jt_totale)
    
    Jr_totale_err = standErr(Jr_totale)
    Jp_totale_err = standErr(Jp_totale)
    Jt_totale_err = standErr(Jt_totale)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_gradNe = hist.Jr_gradNes*1000.
    Jp_gradNe = hist.Jp_gradNes*1000.
    Jt_gradNe = hist.Jt_gradNes*1000.
    
    Jr_gradNe_mean = np.nanmean(Jr_gradNe)
    Jp_gradNe_mean = np.nanmean(Jp_gradNe)
    Jt_gradNe_mean = np.nanmean(Jt_gradNe)
    
    Jr_gradNe_median = np.nanmedian(Jr_gradNe)
    Jp_gradNe_median = np.nanmedian(Jp_gradNe)
    Jt_gradNe_median = np.nanmedian(Jt_gradNe)
    
    Jr_gradNe_err = standErr(Jr_gradNe)
    Jp_gradNe_err = standErr(Jp_gradNe)
    Jt_gradNe_err = standErr(Jt_gradNe)
    
    #Current results scaled to nA/m^2 from their stored values of micro A/m^2
    Jr_divTe = hist.Jr_divTes*1000.
    Jp_divTe = hist.Jp_divTes*1000.
    Jt_divTe = hist.Jt_divTes*1000.
    
    Jr_divTe_mean = np.nanmean(Jr_divTe)
    Jp_divTe_mean = np.nanmean(Jp_divTe)
    Jt_divTe_mean = np.nanmean(Jt_divTe)
    
    Jr_divTe_median = np.nanmedian(Jr_divTe)
    Jp_divTe_median = np.nanmedian(Jp_divTe)
    Jt_divTe_median = np.nanmedian(Jt_divTe)
    
    Jr_divTe_err = standErr(Jr_divTe)
    Jp_divTe_err = standErr(Jp_divTe)
    Jt_divTe_err = standErr(Jt_divTe)
    
    #Prints out data table:
    print("Number of Events:", len(mpc.MP_ID))
    
    print(" ")
    print("Jr_curl: %.2f (%.2f) " %(Jr_curl_mean,Jr_curl_median)+ u"\u00B1 %.2f"%Jr_curl_err)
    print("Jr_perp: %.2f (%.2f) " %(Jr_perp_mean,Jr_perp_median)+ u"\u00B1 %.2f"%Jr_perp_err)
    print("Jr_para: %.2f (%.2f) " %(Jr_para_mean,Jr_para_median)+ u"\u00B1 %.2f"%Jr_para_err)
    
    print(" ")
    print("Jr_totali: %.2f (%.2f) " %(Jr_totali_mean,Jr_totali_median)+ u"\u00B1 %.2f"%Jr_totali_err)
    print("Jr_gradNi: %.2f (%.2f) " %(Jr_gradNi_mean,Jr_gradNi_median)+ u"\u00B1 %.2f"%Jr_gradNi_err)
    print("Jr_divTi: %.2f (%.2f) " %(Jr_divTi_mean,Jr_divTi_median)+ u"\u00B1 %.2f"%Jr_divTi_err)
    
    print(" ")
    print("Jr_totale: %.2f (%.2f) " %(Jr_totale_mean,Jr_totale_median)+ u"\u00B1 %.2f"%Jr_totale_err)
    print("Jr_gradNe: %.2f (%.2f) " %(Jr_gradNe_mean,Jr_gradNe_median)+ u"\u00B1 %.2f"%Jr_gradNe_err)
    print("Jr_divTe: %.2f (%.2f) " %(Jr_divTe_mean,Jr_divTe_median)+ u"\u00B1 %.2f"%Jr_divTe_err)
    print(" ")
    
    print("Jp_curl: %.2f (%.2f) " %(Jp_curl_mean,Jp_curl_median)+ u"\u00B1 %.2f"%Jp_curl_err)
    print("Jp_perp: %.2f (%.2f) " %(Jp_perp_mean,Jp_perp_median)+ u"\u00B1 %.2f"%Jp_perp_err)
    print("Jp_para: %.2f (%.2f) " %(Jp_para_mean,Jp_para_median)+ u"\u00B1 %.2f"%Jp_para_err)
    
    print(" ")
    print("Jp_totali: %.2f (%.2f) " %(Jp_totali_mean,Jp_totali_median)+ u"\u00B1 %.2f"%Jp_totali_err)
    print("Jp_gradNi: %.2f (%.2f) " %(Jp_gradNi_mean,Jp_gradNi_median)+ u"\u00B1 %.2f"%Jp_gradNi_err)
    print("Jp_divTi: %.2f (%.2f) " %(Jp_divTi_mean,Jp_divTi_median)+ u"\u00B1 %.2f"%Jp_divTi_err)
    
    print(" ")
    print("Jp_totale: %.2f (%.2f) " %(Jp_totale_mean,Jp_totale_median)+ u"\u00B1 %.2f"%Jp_totale_err)
    print("Jp_gradNe: %.2f (%.2f) " %(Jp_gradNe_mean,Jp_gradNe_median)+ u"\u00B1 %.2f"%Jp_gradNe_err)
    print("Jp_divTe: %.2f (%.2f) " %(Jp_divTe_mean,Jp_divTe_median)+ u"\u00B1 %.2f"%Jp_divTe_err)
    print(" ")
    
    print("Jt_curl: %.2f (%.2f) " %(Jt_curl_mean,Jt_curl_median)+ u"\u00B1 %.2f"%Jt_curl_err)
    print("Jt_perp: %.2f (%.2f) " %(Jt_perp_mean,Jt_perp_median)+ u"\u00B1 %.2f"%Jt_perp_err)
    print("Jt_para: %.2f (%.2f) " %(Jt_para_mean,Jt_para_median)+ u"\u00B1 %.2f"%Jt_para_err)
    
    print(" ")
    print("Jt_totali: %.2f (%.2f) " %(Jt_totali_mean,Jt_totali_median)+ u"\u00B1 %.2f"%Jt_totali_err)
    print("Jt_gradNi: %.2f (%.2f) " %(Jt_gradNi_mean,Jt_gradNi_median)+ u"\u00B1 %.2f"%Jt_gradNi_err)
    print("Jt_divTi: %.2f (%.2f) " %(Jt_divTi_mean,Jt_divTi_median)+ u"\u00B1 %.2f"%Jt_divTi_err)
    
    print(" ")
    print("Jt_totale: %.2f (%.2f) " %(Jt_totale_mean,Jt_totale_median)+ u"\u00B1 %.2f"%Jt_totale_err)
    print("Jt_gradNe: %.2f (%.2f) " %(Jt_gradNe_mean,Jt_gradNe_median)+ u"\u00B1 %.2f"%Jt_gradNe_err)
    print("Jt_divTe: %.2f (%.2f) " %(Jt_divTe_mean,Jt_divTe_median)+ u"\u00B1 %.2f"%Jt_divTe_err)
    print(" ")
    
    return

def InitializeData():
    
    #Averaged data files: 
    quiet_average_data = "quiet_mpcData.txt"
    edr_average_data = "edr_mpcData.txt"
    
    #Histogram data files:
    quiet_hist_data = "quiet_histData.txt"
    edr_hist_data = "edr_histData.txt"
    
    #Initializes EDR and Quiet classes for the averaged data:
    quiet = magnetopauseCrossingData(225)
    edr = magnetopauseCrossingData(23)
    
    #Initializes EDR and MPC (quiet) classes for the hist data:
    quiet_hist = magnetopauseCrossingData(225)
    edr_hist = magnetopauseCrossingData(23)
    
    #Reads in data from averaged crossings:
    averaged_data_reader(quiet,quiet_average_data)
    averaged_data_reader(edr,edr_average_data)
        
    #Reads in hist data from crossings:
    hist_reader(quiet_hist,quiet_hist_data)    
    hist_reader(edr_hist,edr_hist_data)    
    
    #Plots Current and Velocity histograms
    createHistograms(edr_hist,quiet_hist)
    
    #Creates Table entry for Quiet Events: 
    print(" ")
    print("----------------------------")
    print("Quiet Event Data:")
    data_table(quiet,quiet_hist)
    
    #Creates Table entry for EDR Events: 
    print(" ")
    print("----------------------------")
    print("EDR Event Data:")
    data_table(edr,edr_hist)
    
    return

InitializeData()





